POP Steering parameterization

Eddy diffusivity implementation by Bates, Tulloch, Marshall, and Ferrari

POP Steering Parameterization

CAM Implementation of Mesoscale eddy diffusivity based on mixing theory

Here are my notes on the development of a parameterization for eddy diffusivity. The scheme is based on the paper by Michael Bates, Ross Tulloch, John Marshall, and Raffaele Ferrari. I'm working with Gokhan Danabasoglu and Matt Long here at NCAR and John Marshall (co-author) at MIT.

The general approach of this effort will be to follow Gokhan's implementation notes, isolating each of the terms contributing to the eddy diffusivity given by equation 6 in Bates et al. and comparing the CESM implementation of those terms with Bates, reanalysis, or other observations.

Gokhan's Implementation Notes on Diffusivity with Steering Level Suppression

12 November 2014

Starting with $$K=u_{rms}∗L_{mix}$$

the general form for diffusivity, $K$, is given by equation (6) of Bates et al. (2014). Namely,

$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

Here,

  • $u_{rms}$ is the root-mean-square (rms) eddy velocity;

  • $L_{mix}$ is the mixing length;

  • $L_{eddy}$ is the eddy diameter (depth independent);

  • $u_{mean}$ is the mean zonal velocity (resolved);

  • $c$ is the zonal eddy phase speed (depth independent);

  • $\Gamma = 0.35$;

  • $b1 \sim 4$.

Eddy Length Scales:

Rossby deformation radius = $L_r = {c_r \over |f|}$

Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$

Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

$c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;

$f$ is the Coriolis parameter;

$\beta$ is the latitudinal variation of the Coriolis parameter; and

$\sigma_{vi}$ is the Eady growth rate given by

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

So, any one these length scales could be used as an eddy length scale. An alternative is

$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$

Eddy Velocity:

$u_{rms} = alpha*\sigma_{vi}*L_r$

where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.

Zonal phase speed:

$c = - \beta * L_r^2$

References:

Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44, 1523-1540, doi: 10.1175/JPO-D-13-0130.1.

Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.

Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface altimetric observations in terms of planetary waves and geostropic turbulence. J. Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.

Questions:

  1. In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?

    A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only

  2. In the 2D implementation, vertical profile will be specified by $N2(z)$?

    A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise

  3. Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?

    A: No, $u_{rms}$ is depth independent

  4. alpha = ?

    A: trial and error

  5. Cancellation of $f$’s in $u_{rms}$ calculation?

  6. Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).

    A: What is plotted in Bates is (U-c)

Parameterizing the Eddy Length Scale $\quad \color{red}{L_{eddy}}$

$$K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

The three length scales under consideration:

  1. Rossby Deformation Radius $\quad\quad L_r = {c_r \over |f|}$
  2. Rossby Equatorial Radius $ \quad\quad= L_{req} = {\sqrt{c_r \over 2\beta}}$
  3. Rhine's Scale $ \quad\quad L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

The Rossby Deformation Radius $L_r$ depends on the the Baroclinic wave speed $c_r$ and the Coriolis force $f$. The CESM value of the first Baroclinic wave speed is derived as per eq 2.2 in Chelton 1998.

First Baroclinic Wave Speed $c_r$ Chelton vs CESM

In [39]:
Expand Code

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$L_r = {c_r \over |f|}$

$L_{req} = {\sqrt{c_r \over 2\beta}}$

$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$

$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$

In [40]:
Expand Code
In [41]:
Expand Code

Parameterizing Zonal Eddy Phase Speed $\color{red}{c}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} $

$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator

$c = - \beta * L_{eddy}^2$

In [42]:
Expand Code

Here is the new c_eddy field ( Using $c = max(- \beta * L_r^2,-20)$ )

In [43]:
Expand Code

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

In [44]:
Expand Code

Parameterizing the Zonal Velocity $\color{red}{u_{mean}}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0))} $

In [45]:
Expand Code

${(U-c)}$ (using $c = - \beta * L_r^2$)

In [46]:
Expand Code

${(U-c)}$ ( Using $c = max(- \beta * L_r^2,-20)$ )

In [47]:
Expand Code
In [48]:
Expand Code

$(U-c)^2$ CESM vs Bates et al

$\quad\quad$ First Calculating $(U-c)^2$ ($c = - \beta * L_r^2$)

In [49]:
Expand Code

$\quad\quad (U-c)^2$ using the new c_eddy field ($c = max(- \beta * L_r^2,-20)$ )

In [50]:
Expand Code
In [51]:
Expand Code

Parameterizing the root-mean-square eddy velocity $\color{red}{u_{rms}}$ and the Eady Growth Rate $\color{red}{\sigma_{vi}}$

$K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $

$u_{rms} = alpha*\color{red}{\sigma_{vi}}*L_{eddy}$

The original derivation of $\sigma_{vi}$:

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

An alternate derivation of $\sigma_{vi}$ is:

$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$

$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $

After hydrostatic and geostrophic approximations

$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

so

$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\therefore$

$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$

$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$

$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$

$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$

Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$

$\quad\quad$ so $\cdots$ this is not $R_i$

Implementation notes

Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$

Denominator :

$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $

$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 / DYT(i,j)^2 \\ \end{align} $

$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$

Notes:

1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
In [52]:
Expand Code

$u_{rms} = alpha * \sigma_{new} * l_{eddy} (alpha=4) \quad$ CESM vs Bates

$u_{rms}$ is limited to 5cm/s $ \quad max(u_{rms},5.)$

In [53]:
Expand Code
In [54]:
Expand Code

$u_{rms}^2 \quad$ CESM vs Bates

In [55]:
Expand Code
In [56]:
Expand Code

${(U-c)}^2 \over u_{rms}^2 \quad$ CESM vs Bates

Using $c = max(- \beta * L_r^2,-20)$ and $u_{rms} = max(u_{rms},5.)$

In [57]:
Expand Code

The CESM values are high compared to Bates.

In [58]:
Expand Code

Suppression factor $= {1 \over (1 + b1 * |\bar u - c|^2 /u_{rms (z=0)}^2 )}$

(Adjustment to U_rms helps the suppression everywhere except N/S of 50 degrees.)

In [59]:
Expand Code
In [60]:
Expand Code

Parameterizing the Mixing term $\color{red}{L_{mix}}$

$\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression$

$\Gamma = 0.35$

$b1 = 4$

$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$

$\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

LMIX2

1.new c fixes values around equator
2.latest urms average using m^2/n implementation of sigma
In [61]:
Expand Code

Parameterizing Eddy diffusivity $\color{red}{K}$

$\color{red}{K}=u_{rms}∗L_{mix}$

Shown are

$K$

$K_{lim}$ using a supression field with $u_{rms}^2$ limited to 100.

They both look very similar but the $K_{lim}$ does add more suppression (albeit small) outside of 20 lat. The high maximum value of this field though still causes the model to blow up.

In [62]:
Expand Code
In [63]:
Expand Code
In [64]: 
knosupp = nc.variables['KAPPA_NO_SUPP'][0,:,:]
knosupp_nearest = pyresample.kd_tree.resample_nearest(orig_def,knosupp*1.e-4, 
        targ_def0, radius_of_influence=500000, fill_value=None)

clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(targ_def0,knosupp_nearest,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ With a suppression scaling factor of 1 (max = %.1f/min = %.1f) "%(knosupp_nearest.max(),knosupp_nearest.min()))
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,10e3])
pylab.xticks([0,5e3,10e3]);

Steering Level status

1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale ok
    a. Use alternate calculation otherwise has problems at equator and
    b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c) 
    a. Too high 10S:10N (applying an upper limit of 20cm/s helps)
    b. Missing westward velocities in ACC which appear 
        in Hughes Data used in this parameterization
5. (U-c) Is close if c is limited to 20cm/s
6. U_rms  (alpha*sigma*l_eddy)  OK but low values of U_rms is driving 
   errors in the suppression term.     
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression OK
9. K is now in the ball park but max values are too great.
In [65]:
Expand Code

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters

In [66]:
Expand Code
In [67]:
Expand Code
In [68]:
Expand Code
Out[68]:

Posted by John Truesdale